function out = getAVarTheta11Step3(namesTheta12,theta11Step3,AVarTheta1Step1,AVarTheta2Step2,AVarTheta12Step3,...
    setupStep1,setupStep3,epsValue,theta12OptimalOn)

% The standard errors for the third step WITHOUT correction for estimation uncertainty in theta12
VTheta11Step3  = getAVarTheta1andFactors(theta11Step3,setupStep3,epsValue);

% Computing PSI for theta12 by adding elements of sigma into the set of theta1
if setupStep3.AVarTheta1On == 1
    numTheta12  = size(namesTheta12,1);
    theta1Step3 = theta11Step3;
    for i=1:numTheta12
        theta1Step3.(namesTheta12{i}) = getfield(setupStep3.calibrateTheta1,namesTheta12{i});
    end
    % Note we deliberately use setupStep1
    setupStep1.AVarTheta1On      = 1;
    AVarTheta1Step3              = getAVarTheta1andFactors(theta1Step3,setupStep1,epsValue);
    posTheta12Step3              = zeros(1,numTheta12);
    for i=1:numTheta12
        posTheta12Step3(1,i) = find(strcmp(AVarTheta1Step3.names,namesTheta12(i)));
    end
    
    % Estimation of K
    numTheta11 = size(fieldnames(theta11Step3),1);
    J = zeros(numTheta11,numTheta12);
    T = size(VTheta11Step3.PSI,2);
    nyIndex = zeros(T,1);
    for t=1:T
        nyIndex(t,1) = sum(~isnan(setupStep3.yields(t,:)));
        for j=1:nyIndex(t,1)
            J = J + VTheta11Step3.PSI(:,t,j)*AVarTheta1Step3.PSI(posTheta12Step3,t,j)';
        end
    end
    J = J/sum(nyIndex,1);
    K = - VTheta11Step3.invA_Theta1 * J;
    
    % Estimation of U = Cov(theta12Step1,theta11Step3)
    posTheta11Step3 = zeros(1,numTheta11);
    for i=1:numTheta11
        posTheta11Step3(1,i) = find(strcmp(VTheta11Step3.names,VTheta11Step3.names(i)));
    end
    omega21Step1  = AVarTheta1Step1.invA_Theta1(posTheta12Step3,posTheta11Step3);
    omega22Step1  = AVarTheta1Step1.invA_Theta1(posTheta12Step3,posTheta12Step3);
    URobust       = zeros(numTheta12,numTheta11);
    UHomo         = zeros(numTheta12,numTheta11);
    UHomoPooled   = zeros(numTheta12,numTheta11);
    if theta12OptimalOn == 1
        idx           = 0;
        for t=1+setupStep3.wT:T-setupStep3.wT
            for j=1:nyIndex(t,1)-setupStep3.wD
                idx = idx + 1;
                for kT=-setupStep3.wT:setupStep3.wT
                    for kD=0:setupStep3.wD
                        if kD == 0 % Because then only one term, i.e. no symmetry
                            tmpRobust = (1-abs(kT)/(1+setupStep3.wT))*(1-abs(kD)/(1+setupStep3.wD))...
                                *AVarTheta1Step1.residuals(j,t)*VTheta11Step3.residuals(j+kD,t+kT)...
                                *(omega21Step1*AVarTheta1Step1.PSI(posTheta11Step3,t,j) + ...
                                 omega22Step1*AVarTheta1Step1.PSI(posTheta12Step3,t,j))...
                                 *VTheta11Step3.PSI(:,t+kT,j+kD)'*VTheta11Step3.invA_Theta1;
                        else
                            tmpRobust = (1-abs(kT)/(1+setupStep3.wT))*(1-abs(kD)/(1+setupStep3.wD))...
                                *AVarTheta1Step1.residuals(j,t)*VTheta11Step3.residuals(j+kD,t+kT)...
                                *( (omega21Step1*AVarTheta1Step1.PSI(posTheta11Step3,t,j)    + omega22Step1*AVarTheta1Step1.PSI(posTheta12Step3,t,j))*VTheta11Step3.PSI(:,t+kT,j+kD)'*VTheta11Step3.invA_Theta1 ...
                                +  (omega21Step1*AVarTheta1Step1.PSI(posTheta11Step3,t,j+kD) + omega22Step1*AVarTheta1Step1.PSI(posTheta12Step3,t,j+kD))*VTheta11Step3.PSI(:,t+kT,j)'*VTheta11Step3.invA_Theta1);
                        end
                        if ~isnan(tmpRobust)
                            URobust = URobust + tmpRobust;
                        end
                        
                        if kT == 0 && kD == 0
                            tmpHomo =  VTheta11Step3.RvHat(t,1)...
                                *(omega21Step1*AVarTheta1Step1.PSI(posTheta11Step3,t,j) + omega22Step1*AVarTheta1Step1.PSI(posTheta12Step3,t,j))...
                                *VTheta11Step3.PSI(:,t+kT,j+kD)'*VTheta11Step3.invA_Theta1;
                            tmpHomoPooled = VTheta11Step3.RvHatPooled ...
                                *(omega21Step1*AVarTheta1Step1.PSI(posTheta11Step3,t,j) + omega22Step1*AVarTheta1Step1.PSI(posTheta12Step3,t,j))...
                                *VTheta11Step3.PSI(:,t+kT,j+kD)'*VTheta11Step3.invA_Theta1;
                        else
                            tmpHomo = (1-abs(kT)/(1+setupStep3.wT))*(1-abs(kD)/(1+setupStep3.wD))...
                                *AVarTheta1Step1.residuals(j,t)*VTheta11Step3.residuals(j+kD,t+kT)...
                                *( (omega21Step1*AVarTheta1Step1.PSI(posTheta11Step3,t,j) + omega22Step1*AVarTheta1Step1.PSI(posTheta12Step3,t,j))*VTheta11Step3.PSI(:,t+kT,j+kD)'*VTheta11Step3.invA_Theta1 ...
                                +(omega21Step1*AVarTheta1Step1.PSI(posTheta11Step3,t,j+kD) + omega22Step1*AVarTheta1Step1.PSI(posTheta12Step3,t,j+kD))*VTheta11Step3.PSI(:,t+kT,j)'*VTheta11Step3.invA_Theta1);
                            tmpHomoPooled = tmpHomo;
                        end
                        if ~isnan(tmpHomo)
                            UHomo       = UHomo + tmpHomo;
                            UHomoPooled = UHomoPooled + tmpHomoPooled;
                        end
                    end
                end
            end
        end
        URobust = URobust/idx^2;
        UHomo   = UHomo/idx^2;
    end
    
    % Position of theta12 at step 2
    posTheta12Step2          = zeros(1,numTheta12);
    for i=1:numTheta12
        posTheta12Step2(1,i) = find(strcmp(AVarTheta2Step2.names,namesTheta12(i)));
    end
    
    AVarTheta11Step3Robust     = VTheta11Step3.VarRobust     + K*AVarTheta2Step2.VarRobust(posTheta12Step2,posTheta12Step2)*K' ...
        + K*AVarTheta12Step3.lambdaRobust*URobust ...
        + URobust'*AVarTheta12Step3.lambdaRobust'*K';
    AVarTheta11Step3Homo       = VTheta11Step3.VarHomo       + K*AVarTheta2Step2.VarRobust(posTheta12Step2,posTheta12Step2)*K' ...
        + K*AVarTheta12Step3.lambdaHomo*UHomo ...
        + UHomo'*AVarTheta12Step3.lambdaHomo'*K';
    AVarTheta11Step3HomoPooled = VTheta11Step3.VarHomoPooled + K*AVarTheta2Step2.VarRobust(posTheta12Step2,posTheta12Step2)*K' ...
        + K*AVarTheta12Step3.lambdaHomoPooled*UHomoPooled ...
        + UHomoPooled'*AVarTheta12Step3.lambdaHomoPooled'*K';
else
    numTheta11                 = size(fieldnames(theta11Step3),1);
    AVarTheta11Step3Robust     = NaN(numTheta11,numTheta11);
    AVarTheta11Step3Homo       = NaN(numTheta11,numTheta11);
    AVarTheta11Step3HomoPooled = NaN(numTheta11,numTheta11);
end

% Saving output
out.VarRobust        = AVarTheta11Step3Robust;
out.VarHomo          = AVarTheta11Step3Homo;
out.VarHomoPooled    = AVarTheta11Step3HomoPooled;
out.VarFactorsRobust = VTheta11Step3.VarFactorsRobust;
out.VarFactorsHomo   = VTheta11Step3.VarFactorsHomo;
out.autoCov10Factors = VTheta11Step3.autoCov10Factors;
out.autoCov01Factors = VTheta11Step3.autoCov01Factors;
out.autoCov1Residuals= VTheta11Step3.autoCov1Residuals;
out.varResiduals     = VTheta11Step3.varResiduals;
out.RvHat            = VTheta11Step3.RvHat;
out.RvHat_se         = VTheta11Step3.RvHat_se;
out.RvHatPooled      = VTheta11Step3.RvHatPooled;
out.RvHatPooled_se   = VTheta11Step3.RvHatPooled_se;
out.residuals        = VTheta11Step3.residuals;
out.names            = setupStep3.selectTheta1;
out.PSI              = VTheta11Step3.PSI;
out.A_Theta11        = VTheta11Step3.A_Theta1;

end